
#Clear R
rm(list=ls())

#load necessary function
zero1 <- function(x, minx=NA, maxx=NA){
  res <- NA
  if(is.na(minx)) res <- (x - min(x,na.rm=T))/(max(x,na.rm=T) -min(x,na.rm=T))
  if(!is.na(minx)) res <- (x - minx)/(maxx -minx)
  res
}


#function to install packages if they don't exist
ipak <- function(pkg){
  new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
  if (length(new.pkg)) 
    install.packages(new.pkg, dependencies = TRUE)
  sapply(pkg, require, character.only = TRUE)
}

# usage
packages <- c("ggplot2", "psych","interplot","stargazer","stringr","tidyr","lavaan","dplyr","data.table")
ipak(packages)

load("wvs_sub.RData")

##########################
# Descriptives reported in the paper
#########################

#lr placement
summary(wvs_sub$wvs_lrplacement)
sd(wvs_sub$wvs_lrplacement, na.rm=T)

#cultural ideology
data_social <- (data.frame(wvs_sub$purity, wvs_sub$gay1, wvs_sub$gay2, wvs_sub$abortion, wvs_sub$euthanasia,wvs_sub$immigration))
cronbach(data_social)

summary(wvs_sub$wvs_moral)
sd(wvs_sub$wvs_moral, na.rm=T)

#Economic ideology
data_econ <- (data.frame(wvs_sub$redistribution1, wvs_sub$redistribution2, wvs_sub$govrespon))
cronbach(data_econ)

summary(wvs_sub$wvs_economicattitudes)
sd(wvs_sub$wvs_economicattitudes, na.rm=T)

#correlation between social and economic attitudes: 0.08
cor(data<-(data.frame(wvs_sub$wvs_economicattitudes, wvs_sub$wvs_moral)), use="complete.obs", method="pearson") 
cor(data<-(data.frame(wvs_sub$wvs_moral, wvs_sub$wvs_lrplacement)), use="complete.obs", method="pearson") 
cor(data<-(data.frame(wvs_sub$wvs_economicattitudes, wvs_sub$wvs_lrplacement)), use="complete.obs", method="pearson") 


#########################
#
#Create Figure 5
#
#########################

# BFI models
#moral attitudes
t1 <- lm(wvs_moral~BFI_open+BFI_con+BFI_ext+BFI_agre+BFI_neu+gender+age+age2+as.factor(education)+as.factor(income_categories),data=wvs_sub)
#economic attitudes
t2 <- lm(wvs_economicattitudes~BFI_open+BFI_con+BFI_ext+BFI_agre+BFI_neu+gender+age+age2+as.factor(education)+as.factor(income_categories),data=wvs_sub)
#Left right placement
t3 <- lm(wvs_lrplacement~BFI_open+BFI_con+BFI_ext+BFI_agre+BFI_neu+gender+age+age2+as.factor(education)+as.factor(income_categories),data=wvs_sub)

# Mini IPIP
#moral attitudes (gay rights)
mi1 <- lm(wvs_moral~Mini_IPIP_open+Mini_IPIP_con+Mini_IPIP_ext+Mini_IPIP_agre+Mini_IPIP_neu+gender+age+age2+as.factor(education)+as.factor(income_categories), data=wvs_sub)
#economic attitudes
mi2 <- lm(wvs_economicattitudes~Mini_IPIP_open+Mini_IPIP_con+Mini_IPIP_ext+Mini_IPIP_agre+Mini_IPIP_neu+gender+age+age2+as.factor(education)+as.factor(income_categories), data=wvs_sub)
#Left right placement
mi3 <- lm(wvs_lrplacement~Mini_IPIP_open+Mini_IPIP_con+Mini_IPIP_ext+Mini_IPIP_agre+Mini_IPIP_neu+gender+age+age2+as.factor(education)+as.factor(income_categories), data=wvs_sub)

# Full IPIP
#moral attitudes
i1 <- lm(wvs_moral~IPIP_open+IPIP_con+IPIP_ext+IPIP_agre+IPIP_neu+gender+age+age2+as.factor(education)+as.factor(income_categories), data=wvs_sub)
#economic attitudes
i2 <- lm(wvs_economicattitudes~IPIP_open+IPIP_con+IPIP_ext+IPIP_agre+IPIP_neu+gender+age+age2+as.factor(education)+as.factor(income_categories), data=wvs_sub)
#Left right placement
i3 <- lm(wvs_lrplacement~IPIP_open+IPIP_con+IPIP_ext+IPIP_agre+IPIP_neu+gender+age+age2+as.factor(education)+as.factor(income_categories), data=wvs_sub)
summary(i3)


#Cultural ideology
moralattitudes <- data.frame(rbind(summary(t1)$coefficients[c(2:6),1:2],summary(i1)$coefficients[c(2:6),1:2],summary(mi1)$coefficients[c(2:6),1:2]))
moralattitudes$trait <- c("Openness","Conscientiousness","Extraversion","Agreeableness","Neuroticism")
moralattitudes$scale <- c(rep("BFI",5),rep("IPIP",5),rep("Mini-IPIP",5))
colnames(moralattitudes)[2]="se"

# economic attitudes
economicattitudes <- data.frame(rbind(summary(t2)$coefficients[c(2:6),1:2],summary(i2)$coefficients[c(2:6),1:2],summary(mi2)$coefficients[c(2:6),1:2]))
economicattitudes$trait <- c("Openness","Conscientiousness","Extraversion","Agreeableness","Neuroticism")
economicattitudes$scale <- c(rep("BFI",5),rep("IPIP",5),rep("Mini-IPIP",5))
colnames(economicattitudes)[2]="se"

# left-right behavior
lrbehavior <- data.frame(rbind(summary(t3)$coefficients[c(2:6),1:2],summary(i3)$coefficients[c(2:6),1:2],summary(mi3)$coefficients[c(2:6),1:2]))
lrbehavior$trait <- c("Openness","Conscientiousness","Extraversion","Agreeableness","Neuroticism")
lrbehavior$scale <- c(rep("BFI",5),rep("IPIP",5),rep("Mini-IPIP",5))
colnames(lrbehavior)[2]="se"

forggplot <- rbind(cbind(lrbehavior,DV='Unidimensional \n Ideology'), cbind(moralattitudes,DV='Cultural \n Ideology'),cbind(economicattitudes,DV='Economic \n Ideology'))

#Create figure
forggplot$scale <- factor(forggplot$scale,levels=c("BFI","Mini-IPIP","IPIP"))
forggplot$trait <- factor(forggplot$trait,levels=(c("Openness","Conscientiousness","Neuroticism","Agreeableness","Extraversion")))

#create plot
ggplot(subset(forggplot),aes(x=scale,y=Estimate))+facet_grid(trait~DV,scales = "free")+geom_pointrange(aes(x=scale,ymin=Estimate-1.96*se,ymax=Estimate+1.96*se))+theme_bw()+theme(axis.text.x = element_text(angle = -45, hjust = 0))+geom_hline(yintercept = 0,lty="dashed")+xlab("Battery")+geom_errorbar(aes(x=scale,ymin=Estimate-1.96*se,ymax=Estimate+1.96*se),width=.2)+geom_point()+theme(text=element_text(size=10))

########################################
#
######### Significance of differences
#
########################################
wvs_sub %>% dplyr::select(contains("BFI"), contains("IPIP"),nomem_encr, gender, age, age2, education,income_categories,wvs_lrplacement,wvs_moral,wvs_economicattitudes) -> fomrmelting

melt(setDT(na.omit(fomrmelting)), measure = list(c(1,6,13),c(2,7,14),c(3,8,12),c(4,9,15),c(5,10,11)),value.name = c("Extraversion","Agreeableness","Conscientiousness","Neuroticism","Openness")) -> melted12
melted12$variable <- factor(melted12$variable,labels=c("BFI","IPIP","mini-IPIP"))
a <- (lm(wvs_lrplacement~Extraversion*variable+Agreeableness*variable+Conscientiousness*variable+Neuroticism*variable+Openness*variable+gender+age+age2+as.factor(education)+as.factor(income_categories),data=melted12))

b <- (lm(wvs_moral~Extraversion*variable+Agreeableness*variable+Conscientiousness*variable+Neuroticism*variable+Openness*variable+gender+age+age2+as.factor(education)+as.factor(income_categories),data=melted12))

c <- (lm(wvs_economicattitudes~Extraversion*variable+Agreeableness*variable+Conscientiousness*variable+Neuroticism*variable+Openness*variable+gender+age+age2+as.factor(education)+as.factor(income_categories),data=melted12))

#######################
#Create figure 6 from main paper
#######################


source("Study3_Figure6_Openness.R")
op_random <- forfig
source("Study3_Figure6_Agreeableness.R")
ag_random <- forfig
source("Study3_Figure6_Conscientiousness.R")
con_random <- forfig
source("Study3_Figure6_Extraversion.R")
ext_random <- forfig
source("Study3_Figure6_Neuroticism.R")
neur_random <- forfig

forfigure <- rbind(op_random,ag_random,con_random,ext_random,neur_random)
forfigure$trait <- factor(forfigure$trait,levels=(c("Openness","Conscientiousness","Neuroticism","Agreeableness","Extraversion")))
forfigure$variable <- car::recode(forfigure$variable,"'Left.Right'='Unidimensional \\n Ideology';'Moral'='Cultural \\n Ideology';'Economic.Ideology'='Economic \\n Ideology'")
forfigure$variable <- factor(forfigure$variable,levels=levels(forfigure$variable)[c(3,1,2)])

ggplot(forfigure,aes(x=as.factor(Index),y=value))+facet_grid(trait~variable,scales = "free_y")+geom_boxplot()+theme(text=element_text(size=18))+theme_bw()+ylab("Estimate")+xlab("Number of Items in Scale")+geom_hline(yintercept=0, lty = "dashed")

##############################
### ALPHA
##############################
ggplot(subset(forfigure,Index!=1 & variable=="Unidimensional \n Ideology"),aes(y=value,x=X1))+geom_point()+theme_bw()+xlab("Cronbach's Alpha")+geom_smooth(method="lm")+ylab("b Coefficient")+facet_grid(trait~Index)

ggplot(subset(forfigure,Index!=1 & variable=="Cultural \n Ideology"),aes(y=value,x=X1))+geom_point()+theme_bw()+xlab("Cronbach's Alpha")+geom_smooth(method="lm")+ylab("b Coefficient")+facet_grid(trait~Index)

ggplot(subset(forfigure,Index!=1 & variable=="Economic \n Ideology"),aes(y=value,x=X1))+geom_point()+theme_bw()+xlab("Cronbach's Alpha")+geom_smooth(method="lm")+ylab("b Coefficient")+facet_grid(trait~Index)

##############################
#
### factor scores
#
##############################

################
#
#Factor structure Big 5 with fixed variances
#
################

#factor score B5: acceptable fit
#factor score B5: acceptable fit
B5<-'Open = ~ NA*open1+ open2+open3+open4+open5+open6+open7+open8+open9+open10 
Con = ~NA*con1+con2+con3+con4+con5+con6+con7+con8+con9+con10
Ext = ~NA*ext1+ext2+ext3+ext4+ext5+ext6+ext7+ext8+ext9+ext10
Agre= ~NA*agre1+agre2+agre3+agre4+agre5+agre6+agre7+agre8+agre9+agre10
Neu= ~NA*neu1+neu2+neu3+neu4+neu5+neu6+neu7+neu8+neu9+neu10
Open~~1*Open
Con~~1*Con
Ext~~1*Ext
Agre~~1*Agre
Neu~~1*Neu'


#assess fit & take into account categorical nature of the variables
fit<-cfa(B5,data=wvs_sub)

itemscores <- inspect(fit,what="std")$lambda

fscores <-  data.frame(parameterestimates(fit,standardized = T)$std.all[1:50])
                       names(fscores)[1] <- "std.all"
fscores$itemnumber <- rep(1:10,5)
fscores$trait <- c(rep("Openness",10),rep("Conscientiousness",10),rep("Extraversion",10),rep("Agreeableness",10),rep("Neuroticism",10))
library(tidyr)
forfigure %>% separate(itemstogether, into = paste("V", 1:10, sep = "_")) -> both
library(dplyr)
itemsforfactor <- dplyr::select(both,starts_with("V_"),18)

for(i in 1:10){
  itemsforfactor[itemsforfactor==i & itemsforfactor$trait=="Openness"] <- fscores$std.all[fscores$itemnumber==i & fscores$trait=="Openness"]
  itemsforfactor[itemsforfactor==i & itemsforfactor$trait=="Conscientiousness"] <- fscores$std.all[fscores$itemnumber==i & fscores$trait=="Conscientiousness"]
  itemsforfactor[itemsforfactor==i & itemsforfactor$trait=="Extraversion"] <- fscores$std.all[fscores$itemnumber==i & fscores$trait=="Extraversion"]
  itemsforfactor[itemsforfactor==i & itemsforfactor$trait=="Agreeableness"] <- fscores$std.all[fscores$itemnumber==i & fscores$trait=="Agreeableness"]
  itemsforfactor[itemsforfactor==i & itemsforfactor$trait=="Neuroticism"] <- fscores$std.all[fscores$itemnumber==i & fscores$trait=="Neuroticism"]
  
}
itemsforfactor <- mapply(as.numeric,itemsforfactor)
factormeans <- rowMeans(itemsforfactor[,1:10],na.rm=T)
forfactor <- data.frame(both,factormeans)
ggplot(subset(forfactor,variable=="Unidimensional \n Ideology" & Index!=10),aes(y=value,x=factormeans))+geom_point()+theme_bw()+xlab("Mean Factor Score")+geom_smooth(method="lm")+ylab("b Coefficient")+facet_grid(trait~Index)
ggsave("b5factor_lrideology.png",width=11,height=8)

ggplot(subset(forfactor,Index!=10 & variable=="Cultural \n Ideology"),aes(y=value,x=factormeans))+geom_point()+theme_bw()+xlab("Mean Factor Score")+geom_smooth(method="lm")+ylab("b Coefficient")+facet_grid(trait~Index)
ggsave("b5factor_culideology.png",width=11,height=8)

ggplot(subset(forfactor,Index!=10 & variable=="Economic \n Ideology"),aes(y=value,x=factormeans))+geom_point()+theme_bw()+xlab("Mean Factor Score")+geom_smooth(method="lm")+ylab("b Coefficient")+facet_grid(trait~Index)
ggsave("b5factor_econideology.png",width=11,height=8)

